{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sisl\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\n",
    "\\newcommand\\HH{\\mathbf{H}}\n",
    "\\newcommand\\SO{\\mathbf{S}}\n",
    "\\newcommand\\Scat{\\boldsymbol\\Gamma}\n",
    "\\newcommand\\ket[1]{|#1\\rangle}\n",
    "\\newcommand\\bra[1]{\\langle#1|}\n",
    "\\newcommand\\set[1]{\\{#1\\}}\n",
    "\\newcommand\\eig{\\epsilon}\n",
    "$\n",
    "In this example, you will perform transport calculations by projecting scattering states onto molecular orbitals.  \n",
    "\n",
    "We will select Carbon chains as our electrodes and the C$_{60}$ fullerene as our molecule (it has nice symmetries).\n",
    "The molecular projection method is created from a subset molecule space ($\\set M$):\n",
    "\\begin{equation}\n",
    "  \\HH_{\\set M} \\ket{i} = \\eig_i^{\\set M}\\SO_{\\set M} \\ket{i}\n",
    "\\end{equation}\n",
    "with the projectors orthogonalized through the Löwdin transformation\n",
    "\\begin{equation}\n",
    "  \\ket{i'}=\\SO^{1/2}\\ket i\n",
    "\\end{equation}\n",
    "Then the projection of the scattering states will read\n",
    "\\begin{equation}\n",
    "   \\tilde\\Scat = \\sum\\ket{i'}\\bra{i'}\\Scat \\sum \\ket{j'}\\bra{j'}\n",
    "\\end{equation}\n",
    "which completes the basis transformation.\n",
    "\n",
    "These projectors can be performed in TBtrans by defining the $\\set M$ region and defining which scattering matrices should be projected onto.  \n",
    "\n",
    "From the above few equations it should be obvious that to create such a projection it is necessary to define a device region where the scattering matrices are *living* only on the projected region (i.e. the molecule).\n",
    "\n",
    "---\n",
    "\n",
    "Instead of manually setting up the C$_{60}$ molecule, we find the coordinates on this web-page (http://www.nanotube.msu.edu/fullerene/fullerene-isomers.html). The coordinates are stored in the `C60.xyz` file in this directory. Note that when reading this geometry it does not know the orbital distance, so we have to calculate it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C60 = sisl.Geometry.read('C60.xyz')\n",
    "# Calculate the nearest neighbour distance\n",
    "dist = C60.distance(R=5)\n",
    "C60.atoms.atom[0] = sisl.Atom(6, R=dist[0] + 0.01)\n",
    "print(C60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have the molecule, all we need are some electrodes and the connection from the electrodes to the molecule.\n",
    "This electrode is a 2x2 square lattice with transport along the $x$-direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "elec = sisl.Geometry([0] * 3, sisl.Atom(6, R=1.), [1, 1, 10]).tile(2, 0).tile(2, 1)\n",
    "elec.set_nsc(a=3, b=1)\n",
    "H_elec = sisl.Hamiltonian(elec)\n",
    "H_elec.construct(([0.1, 1.1], [0., -1]))\n",
    "H_elec.write('ELEC.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the final device by making the electrode have 2 screening layers, then the C$_{60}$ and finally the right-electrode (equivalently setup to the left part)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "elec_x = elec.tile(3, 0)\n",
    "# Translate to origo\n",
    "C60 = C60.translate(-C60.center(what='xyz'))\n",
    "C60 = C60.translate([-C60.xyz[:, 0].min(), 0, 0])\n",
    "# Do trickery to make sure the coordinates are consecutive along x\n",
    "C60.set_sc([C60.xyz[:, 0].max() + 1., 1., 1.])\n",
    "device = elec_x.append(C60, 0).append(elec_x, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The full device is now created, and we simply need to create the electronic structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = sisl.Hamiltonian(device)\n",
    "H.construct(([0.1, 1.1], [0., -1]))\n",
    "# Correct the C_60 couplings to something different\n",
    "idx_C60 = np.arange(len(elec_x), len(elec_x) + len(C60), dtype=np.int32)\n",
    "for ia in idx_C60:\n",
    "    idx = device.close(ia, R=[0.1, C60.maxR()])[1]\n",
    "    # On-site is already 0, so don't bother doing anything there\n",
    "    # Split idx into C60 couplings and chain couplings\n",
    "    for i in idx:\n",
    "        if i in idx_C60:\n",
    "            H[ia, i] = -1.5\n",
    "        else:\n",
    "            # Coupling to chain\n",
    "            # Since we are only looping atoms in C60\n",
    "            # we have to also set the coupling into C60\n",
    "            # (to assert Hermiticity)\n",
    "            H[ia, i] = 0.1\n",
    "            H[i, ia] = 0.1\n",
    "H.write('DEVICE.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculating projected transmissions is verbose because multiple things are required to be defined. The following flags are required to be specified, carefully read about each of them:\n",
    "- TBT.Atoms.Device\n",
    "- TBT.Projs (name the projection `C60`)\n",
    "- TBT.Proj.<> (create a HOMO and a LUMO projection)\n",
    "- TBT.Projs.T\n",
    "\n",
    "After having completed the `tbtrans` calculation, you will find the additional file `siesta.TBT.Proj.nc` file which contains the transmissions for the projected states.\n",
    "\n",
    "1. Read in the regular `siesta.TBT.nc` file and plot the transmission\n",
    "2. Extend the transmission plot by adding the different projected transmissions (reading the documentation for the `tbtprojncSileTBtrans` is required).  \n",
    "   The naming scheme of projections is a bit complicated, since each projection may be understood as a separate electrode. Therefore, one should request electrodes via 3 names:\n",
    "      1. origin electrode\n",
    "      2. molecule name (in case of more molecules)\n",
    "      3. projection state\n",
    "    So to choose the _electrode_ LUMO from the left electrode, one would request `Left.C60.LUMO`. If you don't want a projection, simply do: `<>.transmission('Left.C60.LUMO', 'Right')`\n",
    "3. Why do the total transmission not match the projected transmissions for the energies they represent? And why may some projections lead to higher transmissions?  \n",
    "   *HINT*: degeneracy of HOMO and LUMO levels\n",
    "4. Adapt the input to include more states in each HOMO, LUMO projection, rerun and check again.  \n",
    "   *HINT*: The $C_{60}$ molecule has 5 HOMO states, 3 LUMO and 3 LUMO$+1$ states.\n",
    "5. **TIME**: Adapt the input to include every HOMO and LUMO state as a separate projection and run again. Compare each transmission."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
